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its  correlation  with  s’.  This  Model  (1)  is  compared  with  the  least  squares 
estimation  Model  (2),  where  the  signal  to  be  estimated  is  directly  part  of 
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— the  model  with  a coefficient  matrix  R.  The  basic  differences  of  these  two 
models  in  the  framework  of  physical  geodesy  are  pointed  out  by  analyzing 
the  validity  of  the  equation 


s'  = R8, 


that  transforms  one  model  into  the  other,  for  different  cases.  For  clarifi- 
cation purposes  least  squares  filtering,  prediction  and  collocation  are  dis- 
cussed separately.  In  filtering  problems  the  coefficient  matrix  R becomes 
the  unit  matrix  and  by  this  the  two  models  become  identical.  For  prediction 
and  collocation  problems  the  relation  s'  = Rs  is  only  fulfilled  in  the  global 
limit  where  s becomes  either  a continuous  function  on  the  earth  or  an  infi- 
nite set  of  spherical  harmonic  coefficients.  ^Applying  Model  (2),  we  see 
that  for  any  finite  dimension  of  s the  operator^equatlons  of  physical  geo- 
desy are  approximated  by  a finite  matrix  relation  whereas  in  Model  (1) 
the  operator  equations  are  applied  in  their  correct  form  on  a continuous, 
approximate  functions  s. 

Both  methods  have  been  applied  in  a numerical  example  where 
spherical  harmonic  coefficients  of  the  geoid  height  were  estimated  from  geoid 
heights  given  in  a global  regular  point  grid  over  Ihe  sphere.  The  results 
show  that  not  only  the  specific  features  of  the  two  least  squares  estimation 
methods  have  to  be  taken  into  consideration  butto  a high  extort  also  the  charac- 
teristics of  the  involved  gravity  quantities  when  a decision  has  to  be  tnsiV> 
which  of  both  methods  should  be  applied. 
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Part  I;  Theory 


1.  Introduction 

The  method  of  least  squares  collocation  was  introduced  into  physical 
geodesy  in  two  principal  presentations  by  T.  Krarup  (1969)  and  H.  Moritz  (1972). 
Meanwhile,  the  method  — in  this  context  denoted  as  "Method  One"  — is  well  es- 
tablished for  different  types  of  approximation  and  adjustment  problems.  Three 
features  for  application  in  physical  geodesy  are  pointed  cut: 


— We  assume  a certain  number  of  measurements  of  a function 
related  to  the  anomalous  potential  of  the  earth's  gravity  field 
is  given,  such  as  gravity  anomalies,  surface  densities,  de- 
flections of  the  vertical  and  so  on.  Then  the  use  of  a harmonic 
and  regular  covariance  expression  in  collocation  will  lead  to  a 
global  approximation  to  this  function  with  minimum  norm.  The 
norm  is  defined  by  the  covariance  function. 

— The  subsequent  application  of  "the  law  of  propagation  of  covar- 
iances" provides  a global  approximation  to  any  desired  function 
related  to  the  earth's  anomalous  potential.  In  addition,  due  to 
the  same  fact,  least  squares  collocation  allows  the  combination 
of  all  types  of  measurements  related  to  the  anomalous  potential. 

— A basic  quality  of  exact  collocation  is  the  reproduction  of  the 
given  data  by  the  approximating  function.  This  property  is  ful- 
filled in  least  squares  collocation  too,  where  the  signal  part  of 
the  measurements  is  reproduced  by  the  approximating  function . 


Unfortunately,  there  are  some  numerical  problems  in  contrast  to  the  theo- 
retical elegance  of  this  method.  Especially  the  necessary  solution  of  a large 
linear  systerirr  with  a dimension  equal  to  the  number  of  observations  can  cause 
severe  difficulties. 

A similar  least  squares  estimation  method,  in  this  context  denoted  as 
"Method  Two",  is  described  in  (P.  Whittle,  1963).  A treatment  of  the  method  in 
geodesy,  for  example,  is  given  by  S.  Lauer  (1971),  and  by  B.  D.  Tapley  and  B. 

E.  Schutz  (1973).  Also  Method  Two  is  called  least  squares  collocation  in  the 
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geodetic  literature,  compare  (H.  Moritz,  1973,  p.  84),  (H.  Moritz  and  K.  P. 
Schwarz,  1973),  and  (K.  P.  Schwarz,  1974). 

A generalisation  in  combining  both  types  of  least  squares  methods  is 
derived  by  H.  Wolf  (1974). 

K.  P.  Schwarz  (1974)  proved  t|at  for  Method  Two  the  solution  of  a linear 
system  with  a dimension  equal  to  the  number  of  observations  can  be  replaced  by 
a solution  of  a system  with  a dimension  equal  to  the  number  of  unknowns.  This 
property  is  of  great  advantage  for  all  overdetermined  problems. 

The  puipose  of  the  present  paper  is  a comparison  of  the  models  under- 
lying these  two  least  squares  estimation  methods  and  an  analysis  of  their  indi- 
vidual features,  since  in  the  geodetic  literature  both  methods  are  not  oleaiiy  discrim- 
inated. lb  addition,  we  expect  from  the  comparison,  an  answer  to  the  question; 
if  the  solution  of  a large  linear  system  in  Method  One,  with  a dimension  equal ' 
to  the  number  of  observations,  can  be  avoided  as  in  Method  Two.  The  considera- 
tions will  be  restricted  to  the  application  of  these  models  to  quantities  related  to 
the  anomalous  potential  of  the  earth. 

hi  a recent  paper  B.  D.  Tapley  (1975)  analyzes  the  two  methods  in  re- 
lation to  the  minimum  variance  principle. 
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A detailed  exposition  of  the  models  and  their  least  squares  solution  is 
given  in  the  above  mentioned  literature.  We  restrict  ourselves  on  the  presenta- 
tion of  the  basic  formulas;  for  Method  One*  equations  (1)  to  (6))  the  corresponding 
equations  for  Method  Two  are  equations  (7)  to  (12). 

The  model  equation  for  Method  One  is 
(1)  j l = Ax  + s'  + n 

j l vector  of  observations,  dim  (n  x 1), 

x unknown  parameters,  dim  (m  x 1), 

A coefficient  matrix,  dim  (n  x m),  that  relates  the 

observations  l to  the  parameters  x, 

s’ random  signal  part  of  l,  dim  (n  x 1),  with  E(s*)  =0 

and  E(s,s,t)  = Ci»*»  , where  E(.)  means  the  statistical 
expectation  and  is  usually  defined  in  physical  geodesy  to 
be  the  integral  over  the  earth,  C*  is  the  covariance 
matrix  for  the  signal  s’. 

n random  signal,  with  E(n)  = 0 and  E(nnT)  = Cu. 

We  also  assume  that  E(s'nT)  = 0. 


The  desired  signal  is  s with  dimension  (k  x 1).  The  signal  s does  not  appear 
explicitly  in  equation  (1),  it  is  linked  to  the  model  by.  a zero  matrix  0.  . We  obtain 


[?]  = B and  [s'® n]  = v , 

t = Ax  + Btv 
where  I . . . identity  matrix. 

The  minimum  length  of  the  vector  v will  be  derived  from  the  minimisation  of 

vT  Q-1v 

where 

,3>  « - &.  S'.  ♦ CJ 

and  C« covariance  vector  of  s 

C«.»  . ...croescovarianoe  between  s and  s'. 

It  is  assumed  that  E(snT)  = 0.  From  the  least  squares  solution  we  obtain 

v = QBt  (BQBt)_1  ( l-  A X) 

and  for  the  parameters  x, 

(4)  x = <At  (BQBTf 1 A)"1  At  (BQBt f1  l. 


where  with  equation  (3) 


c **  + CBB. 


«B'  - ♦ cj  *“d  BQB'  ■ 


Thus,  we  obtain  for  the  components  of  v 


(5)  s - Ci*  (C**  + (t  ■ Ax) 


| 

i 


and 

(6)  s'  + n = (Cjj  + (C**  + CBB)-1  (i-  Ax) 

= l - Ax. 


The  corresponding  derivation  for  Method  Two  starts,  with  the  model 


(7)  l - Ax+Rs  + n 

' ’ u w «i  * ki  a 


Here  the  desired  signal  s appears  explicity  in  the  model  and  is  related  to  the 
observation  vector  l by  the  coefficient  matrix  K with  dim  (n  x k). 

It  has  to  be  emphasized  that  even  though  we  use  in  Model  Two  the  same 
notation  for  x,  s,  and  n these  quantities  may  be  different  from  the  corresponding 
variables  in  Model  Two.  In  fact,  the  purpose  of  the  following  considerations  is  to 
clarify  under  what  circumstances  the  corresponding  variables  x,  s,  and  n are 
identical. 

A rearrangement  of  the  random  quantities  s and  n leads  to 


and  with  B*  = and  v*  = nj  we  obtain  a model  of  the  same  form  as  for 

Method  One 


l = Ax  + B*T  v*. 


We  minimize  v*T  Q*“x  v*  with 


and  obtain  from  the  least  squares  solution 


v*  = Q*  B*t  (B*  Q*  B*t  )"x  ( f - A x) 

and 

(10)  x = (At  (B*Q*B*t)_1  A)'1  At  (B*Q*B*t)-1  l 
where  with  equation  (9) 

Q*B*t  = [C0R]  and  B*Q*B*T  = RC.,RT  + CM. 

Finally,  the  components  of  v*  are 

(11)  s = C„Rt(RCmRt  + C„)‘x  (i-A*) 
and 

(12)  n = Cu(RCh  Rt  + C»0~1(f-Ax) 

With  equations  (1)  to  (6)  we  draw  the  following  conclusions  for  Model  One: 
— The  random  signal  s'  is  of  the  same  physical  nature  as  the  observations  l. 
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-The  desired  vector  s can  be  of  any  dimension,  k,  and  of  any  kind  as  long  as 
it  is  correlated  with  s'.  If  there  exists  no  correlation  between  s and  s', 
which  means  that  there  exists  no  dependence  of  s on  the  observations  A, 
the  covariance  C,',  becomes  zero  and  by  this  also  the  estimate  of  the  sigi- 
nal  s. 


-Although  solution  equation  (6)  for  s'  + n is  trivial  as  it  shows  only  the 
basic  model  (1),  it  nevertheless  reflects  a basic  characteristic  of  all 
collocation  methods,  which  is  the  reproduction  of  the  measurement — in 
our  case  of  the  trendfree  measurement  A - Ax  — at  the  sample  points. 

-If  the  elements  of  the  covariance  matrices  are  built  from  global  covariance 
expressions  fulfilling  Laplace's  equation  and  following  the  'law  of  propa- 
gation" such  as  the  covariance  expressions  derived  by  C.  Tschemlng  and 
R.  Rapp  (1974),  the  solution  equation  (5)  will  converge  for  limit  n->® 
towards  the  linear  operator  equations  in  physical  geodesy,  as  proved  by 
H.  Moritz  (1975). 


From  equations  (7)  and  (8)  for  Model  Two  we  see  that  the  coefficient  matrix 
R has  to  be  given  explicitly. 


r 


3.  Comparison 

A formal  investigation  of  the  equations  for  the  two  models,  (1)  to  (6), 
respectively  (7)  to  (12),  shows  the  following  differences: 


1.  The  least  squares  norm  for  vTv  and  v*fv*  is  obtained  with  respect 
to  different  variance  covariance  matrices  Q,  respectively  Q*. 

2.  The  difference  in  Q and  Q*  causes  the  different  solution  equations 
for  the  trendfree  observation  s'  + n,  equation  (6),  and  for  the 
observation  noise,  n,  equation  (12).  An  equivalent  equation  for 

n can  be  deduced  in  Model  One  by  splitting  equation  (6)  an 
estimation  equation  for  s'  and  one  for  n. 

3.  Due  to  the  same  reason,  i.  e.  the  difference  in  Q,  also  the  solutions 
for  the  parameters  x,  equations  (4)  and  (10),  and  for  the  random 
signal  s,  equations  (5)  and  (11),  are  different  although  the  formulas 
show  formally  a similar  structure. 


The  purpose  of  both  models  is  the  optimal  estimation  of  the  parameters  x and 
of  the  random  signal  s.  Thus,  the  analysis  will  be  concentrated  on  the  solution 
equations  for  x and  s.  We  see  that  for. 


(13)  C*.  = C,.Rt 


and  for 


(14)  C,v  = RC„Rt 


the  solution  equations  (4)  and  (10)  and  also  (5)  and  (11)  become  Identical.  These 
two  expressions  can  be  derived  from  one  single  equation. 


(15)  s’  = RS 


t 


4 
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The  same  relation  is  immediately  obtained  comparing  model  equations  (2)  and  (8). 
The  derivation  of  equations  (13)  and  (14)  from  equation  (15)  is  obvious: 

Equation  (15)  is  multiplied  on  both  sides  with  s and  the  statistical  expectation 
is  taken.  We  obtain: 


(17)  Ets's1)  = RE(sst) 

This  is  a discrete  form  of  the  well  known  WEENER-HOPF  equation,  hi  geodesy  it 
is  the  foundation  of  the  'law  of  propagation  of  covariances*.  Again  the  integral  over 
the  earth  takes  the  place  of  the  statistical  expectation.  The  expectations  E(s'sT) 
and  E (ssT ) become  the  covariances  C,u  and  C„ , and  equation  (17)  becomes 
equation  (13). 

In  the  same  way,  we  take  the  expectation  on  the  square  of  equation  (15) 
and  obtain. 


(18) 


E (s'd*)  = RE(sst)Rt. 


We  replace  the  expectations  by  the  corresponding  covariance  expressions  and  deduce 
equation  (14).  The  comparison  shows  that  the  estimates  for  the  parameters  x and 
the  random  signal  s become  identical  for  both  models,  if  only  equation  (15)  is 
valid.  Therefore  our  further  considerations  will  be  concentrated  on  this  equation. 
We  analyze  if  there  exists  a relation  (15)  with  dimension  k between  two  quantities 
related  to  the  anomalous  potential  of  the  earth.  In  order  to  provide  an  easy  insight 
into  possible  agreements  and  disagreements  the  discussion  is  split  into  filtering, 
prediction  and  collocation. 


3.1  Filtering 

The  purpose  of  least  squares  filtering  is  the  separation  of  the  noise  n from 
the  signal  s in  the  least  squares  sense.  Thereby  we  assume  that  the  covariance 
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function  of  the  signal  and  the  covariance  function  of  signal  plus  noise  are  given. 
Since  filtering  is  usually  carried  out  in  all  observation  points,  we  have  k = n. 
Because  every  signal  s corresponds  to  a certain  observation  i the  coefficient 
matrix  will  be  the  unit  matrix,  R = I.  Thus,  expression  (15)  becomes 


s'  = Is 


or 


(19) 


S'  = 8 


and  identical  estimation  equations  (4)  and  (10),  and  also  (5)  and  (11)  are  ef- 
fective. For  least  squares  filtering  both  models  lead  to  identical  results.. 


3.2  Prediction 


Least  squares  prediction  means  the  optimal  linear  estimation  of  the  ran- 
com  signal  s from  the  observations  i at  k points  which  are  not  identical  with 
the  n observation  points,  hi  physical  geodesy  this  definition  has  to  be  refined. 
We  assume  that  the  observation  and  prediction  points  are  on  the  same  sphere,  <r, 
or  in  spherical  approximation  on  the  earth's  surface  or  cm  the  geoid. 

Following  Krarup  (1969,  p.  16),  we  denote  the  last  two  terms  in  equation 
(5)  with  Knt 


(20) 


in  = (CW  + Cm)-1  ( t - Ax),  i = l,...n. 


The  variable  in , is  only  dependent  on  the  n observation  points  Pt  but  not  on  the 
prediction  points.  The  decision  at  what  and  at  how  many  points  prediction  is  car- 
ried out  depends  on  C„*.  The  covariance  C,*  is  a function  of  the  observation 
points  Pt  and  of  the  finite  or  infinite  number  k of  prediction  points  Q.  The 
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most  general  form  for  the  covariance  vector  C(a<  is  derived  when  each  vector 
element  is  an  analytical  function  of  the  prediction  point  Q,  Qc  o. 


C.i  = C„(Q). 

Then  formula  (5)  becomes  with  (20) 

(21)  s(Q)  = C„(Q)4« 


where  s(Q)  is  a continuous  function  of  the  predicted  signal,  defined  on  o.  This 
case  expresses  the  limit  k-*®  where  s forms  a complete  set. 

Quite  different  is  the  situation  for  Model  Two:  The  corresponding  expres- 
sion for  becomes  for  Method  Two,  (from  (11) ) 

(22)  = (RC„Rr  + C„f 1 (f  - Ax),  J = 1,  ...k. 

i = 1,  ...  n. 

In  contrast  to  equation  (20),  the  expression  £piqj  is  dependent  on  the  observa- 
tion and  prediction  points  since  it  contains  R and  R relates  by  definition,  see  e- 
quation  (7), the  observations  to  the  signal.  Therefore  every  choice  of  k makes 
necessary  a new  computation  of  £Mqj  • 

But  even  if  we  accept  this  restriction  for  £ , the  prediction  cannot  be 
solved  with  Model  Two.  For,  there  does  not  exist  a matrix  R transforming  a 
finite  set  of  random  signals  s into  a set  of  signals  s',  both  quantities  related 
to  the  anomalous  potential  of  the  earth.  Only  for  the  limit  k a solution  is 
theoretically  meaningful.  For  this  limit  die  coefficient  matrix,  R,  degenerates 
towards  Dirac's  functions. 


(23) 


R - ft(^iq) 
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where  4*iq  is  the  spherical  distance  between  the  observation  point  Pt  and  any 
prediction  point  Q.  The  basic  relation  (15)  becomes  for  this  limit 


(24) 


11m  sVi  = -1-  [ 8«p»iq)  S(Q)dCT. 

k -»®  iL 


I 


I 


We  draw  the  conclusions: 

The  least  squares  prediction  problem  can  be  efficiently  treated  by  Model 
One  for  a finite  as  well  as  for  an  infinite  number  of  random  signals.  With  Model 
Two  the  prediction  problem  cannot  be  treated.  The  limit  relation  for  k-«  be- 
tween s'  and  s is  given  in  formula  (24).  Considering  this  expression.  Model 
One  can  be  interpreted  to  be  the  limit  of  Model  Two  for  k , 


3.3  Collocation 


We  distinguish  between  filtering,  prediction,  and  collocation  to  provide 
an  easy  insight  into  die  special  features  of  both  models,  although  in  the  geodetic 
literature  the  term  least  squares  collocation  usually  includes  filtering  and  pre- 
diction. In  least  squares  collocation  gravity  quantities  at  a or  in  E — E 
denotes  the  space  outside  a — are  estimated  from  given  observations  which 
may  be  at  a or  in  E,  Thereby,  the  collocation  quantities  s and  the  observa- 
tions l may  be  gravity  quantities  of  a different  nature  such  as  die  disturbing  po- 
tential, gravity  anomalies,  and  their  first  or  second  derivatives.  For  an  expo- 
sition see  H.  Moritz  (1973),  T.  Krarup  (1969),  and  C.  Tscherning  (1975). 

As  in  least  squares  prediction,'  for  the  more  general  collocation  inethod 
no  finite  relation  (15)  can  be  derived,  now  considering,  as  mentioned  above,  s' 
and  s to  be  gravity  quantities  of  a different  nature.  Only  for  the  extreme  case 
where  s becomes  a continuous  function  oh  o,  the  relation  (15)  between  s'  and  s 
is  satisfied  by  one  of  the  well  known  operator  equations  in  physical  geodesy. 

Example:  The  derivation  of  gravity  anomalies  from  altimeter  data,  as  described 
by  R.  Rapp  (1974), 


Assuming  that  the  altimeter  data  are  synonymous  with  geo  id  un- 
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dulations,  the  observation  vector  1 and  'he  random  signal 
vector  s'  will  consist  of  undulations.  The  random  signal  s 
will  be  a continuous  function  s(Q)  of  the  desired  gravity  ano- 
malies. Then  equation  (15),  l.  e.  the  relation  between  the  un- 
dulations s'  and  the  gravity  anomalies  s is  Stokes’  formula. 


(25) 


with: 


4ttG 


i 


S t (tin q)  s(Q)do 


r. mean  radius  of  the  earth 

G mean  gravity  over  the  earth 

St  (^i) Stokes’  function  for  the  spherical  distance  4>. 


This  means,  in  order  to  fulfill  the  relation  between  s'  and  s 
the  finite  k -dimensional  vector  s becomes  a continuous  func- 
tion defined  on  <r.  Also  in  the  spectral  domain,  in  terms  of  a 
spherical  harmonic  expansion,  only  the  limit  k-®  is  valid. 

The  random  signal  s becomes  the  vector  of  the  spherical  har- 
monic coefficients  of  the  gravity  anomaly  expansion.  We  ob- 
tain for  equation  (15), 

(26)  sVt  = -£*-  £ ■-  - £ ^.(cosv*,)  (s^  cos mXQ  + s£.sinmA„) 

n=a  » = 0 

with: 

PM  (cos  £q)  ....associated  Legendre  polynomial  of  degree  n and  order  m. 


The  equivalence  of  both  models  can  only  be  achieved  for  the  limit  k-» , where 

the  finite  dimensional  matrix,  R,  becomes  one  of  the  well  known  operator  expressions. 


MHMMma 


\ 


For  any  finite  k the  two  models  will  give  different  estimates  for  the  random 
signal  s and  also  for  the  parameters  x. 

Collocation  based  on  Model  One  is.  carried  out  by  inserting  the  continuous 
function  obtained  from  equation  (21)  by  least  squares  prediction  into  an  operator 
equation  that  provides  the  desired  random  signal.  Thereby,  the  continuous  signal 
obtained  from  equation  (21)  which  is  of  the  same  physical  nature  as  the  observa- 
tions is  transformed  into  a finite  or  infinite  number  k of  signals  s. 

Example:  Let  us  again  return  to  the  problem  of  gravity  anomaly  recovery  from 
geoid  undulations. 


From  equation  (21)  an  intermediate  signal  9,  i.e.  a continuous 
undulation  function  defined  on  a is  obtained.  The  covariance 
function  Cpt(Q)  is  in  this  case  the  autocovariance  function 
ciTN)  of  the  geoid  undulations,  hi  order  to  deduce  gravity  ano- 
malies the  signal  9 is  inserted  into  the  boundary  condition  of 
physical  geodesy, 

(27)  8,,  • -G  X, 

where  Sq,  is  the  vector  of  the  gravity  anomalies,  J = 1,  ...  k. 
With  equation  (21)  we  find 

(28)  a,,  = -G  C??(Q)U 

or 


= (-  c*(Q))g€„. 
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We  see  that  in  equation  (28)  the  operator  equation  (27)  is 
only  applied  on  the  covariance  Cff  (Q)  leaving  untouched. 
In  other  words,  the  'law  of  propagation  of  covariances"  is 
applied  on  c£(Q).  The  expression  in  parenthesis  becomes 
the  crosscovariance  between  the  geo  id  undulations  and 
gravity  anomalies,  C^Q).  We  derive  from  equation  (28), 


(29) 


Sqj  = G$n  . 


hi  Model  One  the  only  approximation  is  performed  in  the  prediction  step, 
equation  (21),  where  from  n observations  a function  covering  the  entire  sptere 
o is  estimated.  The  second  step,  where  this  function  is  inserted  into  an  operator 
equation,  is  carried  out  without  any  approximation. 


In  Model  Two  a coefficient  matrix  Rhas  to  be  established.  In  all  practical  situations 
the  dimension k will  be  chosen  to  be  finite  as  shown  for  example  in  (Schwarz,  1974).  In  or- 
der to  solve  a particular  problem,  an  integral  operator  or  its  development  into 
an  infinite  Legendre  polynomial  connecting  the  signal  with  the  observations  has 
to  be  approximated  by  a (n  x k) -dimensional  finite  matrix  R. 

Example:  For  the  application  of  Model  Two  to  our  problem,  the  recovery 

of  gravity  anomalies  from  altimeter  data,  a finite  approximation 
to  equations  (25)  and  (26)  has  to  be  formed.  In  one  case  the 
signal  vector  s will  consist  of  gravity  anomalies  and  the  coeffi- 
cient matrix  R of  n xk  values  of  Stokes'  function.  In  the  other 
case,  the  coefficients  of  a spherical  harmonic  expansion  of  the 
gravity  anomaly  field  are  the  elements  of  the  vector  s and  R 
consists  of  n x k spherical  harmonics. 

In  contrast  to  Model  One,  for  Model  Two  the  operator  equations  have  to  be  approxi- 
mated by  a finite  system. 

Summarizing  the  situation  in  least  squares  collocation,  we  see:  The 
equivalence  between  Model  One  and  Model  Two  can  only  be  achieved  for  the  global 
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limit  k-«.  Then  the  basic  expression  (IS)  becomes  one  of  the  operator  equations 
of  physical  geodesy.  For  any  finite  k collocation  by  Model  One  and  by  Model  Two 
will  result  in  different  estimates  for  s and  for  x.  hi  Model  One,  the  only  approx- 
imation is  introduced  for  the  derivation  of  a function  of  the  type  of  the  observed 
quantities  covering  the  entire  sphere  a,  whereas  in  Model  Two,  one  of  the  inte- 
gral or  spectral  operators  is  approximated  by  a finite  system. 


Part  II:  Numerical  Example 


In  this  second  part  a numerical  example  is  presented  that  allows  us 
to  discuss  some  details  in  the  comparison  of  the  two  least  squares  estimation 
methods.  The  example  treats  the  determination  of  spherical  harmonic  coeffi- 
cients from  geoid  heights.  The  solution  equations  of  the  two  methods  are  un- 
fortunately very  similar  for  this  example  because  of  the  orthogonality  proper- 
ties of  the  spherical  harmonics,  as  will  be  explained  later  in  the  paper.  But 
on  the  other  hand  this  example  has  the  advantage  that  it  can  be  judged  indepen- 
dently of  the  chosen  data.  In  addition,  it  is  to  a certain  extent  complementary 
to  (Schwarz,  1975)  where  for  a similar  problem  regular  least  squares  adjust- 
ment using  no  apriori  covariance  information  is  compared  with  method  (2). 

We  assume  to  have  given  60  geoid  heights  free  of  errors  in  a glo- 
bal spherical  grid  from  latitude  <p  = 60°  to  tfi  * - 60°  and  from  longitude  X = 0° 
to  X * 330°  in  30°  intervals.  The  grid  allows  the  resolution  of  spherical  har- 
monic coefficients  up  to  about  degree  i =*  6 and  order  lm|=  5 . Although  it 
would  theoretically  be  possible  to  determine  from  both  methods,  method  (1) 
as  well  as  method  (2),  an  infinite  number  of  coefficients  — where  the  coef- 
ficients with  l>  6 would  not  have  any  meaning  besides  that  they  would  con- 
tribute to  reproduce  the  original  geoid  heights  — we  solve  only  for  the  coef- 
ficients up  to  t0  * 4 (to  .. . maximum  degree)  and  |m|=  4 . That  means 
that  we  try  to  determine  k = 21  coefficients  (neglecting  degree  zero  and  one) 
given  n = 60  observations. 

The  spherical  harmonic  expansion  of  the  vector,  u,  of  60  geoid 

heights  is, 

(30)  u = ^ ^ 5/.  Yi.  (Pi)  = _T  a , i » 1,  ...,  60; 

* 1 tmo  ■ = -!  a 00  oo X 

with  &U fully  normalized  spherical  harmonic  coefficients 

of  degree  l and  order  m of  the  geoid  height 

a vector  of  these  coefficients  (dim.  «•  1) 


Yj j. (Lilly  normalized  spherical  harmonics 

of  degree  l and  order  m 

r matrix  to  the  harmonics  (dim.  n*°») 


and  with 


where 


Y £.  (Pi)  « (-1)*  [ m)  1 ]*  P£.  (®toLCh)  e1*^ 

60« Kronecker  symbol 

P^(sin^j).  associated  Legendre  functions. 


For  later  needs  the  spherical  harmonic  development  is  splitted  into  a k-dimen- 
tional  part  one  and  a (»  - k + 1)  - dimensional  part  two, 


(31)  u=£a  = R^aj  + S a^ 

»1  b os  ao*  mltkl  B (oo-k  +1)  (co-k  +l)  1 


Method  (I); 

Based  on  the  given  60  geoid  heights  an  approximation  to  the  "true" 
global  geoid  height  function  u(Q)  is  derived  from 


(32)  u(Q)  = Cp„(Q,  PjJC^P,,  Pi)"lu(Pi) 


where  Cj,u  is  the  autocovariance  matrix  of  the  geoid  heights.  The  estimates 
for  the  geoid  height  spherical  harmonic  coefficients  are  obtained  by  inserting 
expression  (32)  into  an  integral  formula, 
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h.  - -jr  J Y«-<Q) 


da 


= | Cuu(Q,Pj)Y£.(Q)da]cttU(P3,Pir1u(PI) 

By  using  for  the  covariance  function  the  expression 
CUU(Q,P)=  £ dj>  s*+a  P£(cos<M 


C Q 

where  dt  = -=V  —ts  * * * 6601(1  hei8ht  degree  variance 

G (*-1) 


square  ratio  of  the  mean  earth 
radius  R,  to  a Bjerhammar 
sphere  with  R^j.  < R, 


we  obtain  for  a^,  , 


(33) 


iu-  (sTTT  s^+1  Y*'  (Pj))  iuu(P3,Pi)1  -(P,) 


and  denoting 


(34) 


fi-  ‘ 


we  find  the  estimation  formula  for  the  geoid  height  coefficients 
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I1  (35) 


Because  of  the  global  character  of  the  approximation  u(Q)  any  geoid  height 
coefficient  may  be  estimated  by  equation  (35)  independently  from  one  another. 
Since  we  solve  only  for  the  lowest  k = 21  coefficients  (zero  and  first  degree 
excluded)  we  use  the  notation  of  equation  (31)  and  rewrite  equation  (35)  to 


11  1 1 Bill 


(36) 


and 

(37) 

Denoting 

and 


we  obtain 


kl  k ■ anal 


a>  = C«u  C'u1  R |j  + C,u  Cuu1  S ap 

kl  k.  « a nkkl  ka  an  a (®-  k + l)  (a>“k  + 1)  1 


Aj 


R 


Aa 


= Cm  Cuu  s 


(38) 


A 


^ Sa  + £2  !h? 


Matrix  relates  the  unknown  low  frequency  spherical  harmonic  coefficients  of 
the  given  geoid  heights  to  the  unknown  spherical  harmonic  coefficients  to  be  esti- 
mated. In  an  ideal  situation  Aj  is  expected  to  be  the  identity  matrix  L Matrix 
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As  relates  the  coefficients  with  l > l0  to  the  coefficients  to  be  estimated. 
Meally  Ao  should  be  the  zero  matrix  O.  In  a practical  situation  like  ours 
where  finite  but  regular  grid  with  observations  defines  a smallest  resolv- 
able wavelength  (or  Nyquist  frequency)  the  matrix  Ap  ji  O generates  the 
alalsing  effect  from  the  coefficients  with  1 > Iq. 


Method  (2) 


Since  we  solve  in  our  case  only  for  k = 21  spherical  harmonic  geo  id 
height  coefficients,  method  (2)  is  based  on  a finite  approximation  to  the  in- 
finite series  expansion  of  equation  (90).  We  express  the  geoid  heights  by  a 
finite  series  of  spherical  harmonics  accepting  a certain  model  error,  c, 


(39)  u = R aj  + f 

ml  akkl  ml 


The  solution  equation  for  the  unknown  parameter  vector  is,  using  equation 
(11) 


(40) 


a j - W (R  C„  ^ 

k l kkkTmkkk  km 


+ S,0'1  « 

am  a 1 


Here  CM  is  the  autocovariance  matrix  of  the  spherical  harmonic  coefficients. 
Because  of  the  orthogonality  of  the  spherical  harmonics  it  becomes  a diagonal 


matrix  with  elements 


1+ i 


and  it  holds  the  relationship 


2L  + 1 


Relation  (41)  based  on  the  orthogonality  properties  of  spherical  harmonics  and 
corresponding  to  equation  (13)  of  part  bne' is  the  reason  that  for  this  specific 
example,  i.  e.  the  determination  of  spherical  harmonic  coefficients  from  a gravi- 
ty function,  the  estimation  equations  of  method  (1)  and  method  (2)  become  very 
similar. 


In  addition,  the  matrix  product  R Q,  RT  is  the  autocovariance  matrix 
for  geoid  height  including  wavelengths  up  t <T  f,,.  We  may  write 

(42)  ci°u  = R C„  Rt 
and 

(43)  Cuu  = Ca,uu  + 


From  equations  (42)  and  (43)  follows  that  relation  (14)  of  part  'one'  does  not 
hold  for  our  example  which  means  that  the  corresponding  estimation  equations 
for  both  methods  are  not  identical. 

We  obtain  for  equation  (40)  using  (41)  and  (42) 

(44)  aj  = C,u  (C^w  + C*,)-1^ 
and  with  (31) 

(45)  J£i  = £»<i  (Ca,«u  + CM)  1 + C*a  (Ca,<m  + C a »)  1Sa?. 


Denoting 


B-s  - C,»  (2a, uu  + Cnn)  S 
the  corresponding  equation  to  (38)  becomes 
(46)  E = 5i  fh  + Sa 

where  again  in  an  ideal  situation  should  be  the  identity  matrix  and  Bp  the 
zero  matrix. 

Although  the  given  geo  id  heights  are  assumed  to  be  free  of  errors 
the  matrix  C„  is  included  in  the  estimation  equation  of  method  (2).  It  may  be 
used  (i)  to  model  partially  the  model  error, €,  in  equation  (39)  and  (ii)  to  sta- 
bilize the  numerical  solution  of  equations  (44)  to  (46).  In  order  to  see  the  in- 
fluence of  the  noise  matrix  CpB  on  the  estimates  two  different  models  were  used: 
model  21  ("£"  for  low  noise  level)  with  £„  = 0.01 1 and  model  2h  ("h"  for  high 
noise  level)  with  CBB  = oa  I with  " 

\ s'*1 

l = Jto+l 


Model  2h  is  corresponding  with  c to  be  interpreted  as  white  noise  with  variance 
equal  to  the  variance  of  the  short  wavelengths  (J  > Iq)  not  contained  in  cj%u . 


Numerical  Results  and  Their  Interpretation 


First,  the  matrices  Ax  for  method  (1)  and  B}  for  method  (2)  have 
been  computed.  They  depend  on  the  spherical  harmonic  coefficients  to  be  esti- 
mated and  on  the  location  of  the  "observations"  but  not  on  their  magnitude.  There- 
fore they  give  an  objective  picture  of  the  transformation  from  a^  to  4j.  For  the 


' 


described  data  grid  the  results  for  matrix  Ax  and  for  matrix  Bx  for  version  21 
(low  noise)  and  2h  (high  noise)  are  given  In  tables  one  to  three. 

The  result  for  matrix  R of  model  21  is  within  1%  a perfect  Identity 
matrix  as  shown  in  Table  2.  Neglecting  the  influence  of  the  second  right  hand 
term  of  equation  (46)  this  means  that  method  (2)  would  produce  perfect  coeffi- 
cient estimates  using  model  21.  Matrix  Bj  (table  three)  using  model  2h  shows 
some  small  off-diagional  terms  and  the  diagonal  terms  slightly  damped  with 
increasing  degree  and  order. 

In  order  to  understand  the  effect  of  Aa  and  Br,  as  well  as  of  the  not 
ideal  matrix  or  Bx  cm  the  estimated  coefficients  we  generated  geoid  height 
spherical  harmonic  coefficients  from  the  GEM  8 set  of  potential  coefficients 
(Wagner  et  als. , 1976)  which  is  complete  up  to  degree  25  with  some  additional 
terms,  and  estimated  the  vector  a^  up  to  = 4 and|m|=  4 using  equations 
(38)  and  (46).  We  split  the  estimated  coefficients  into  three  parts: 

(1)  The  part  produced  by  the  diagonal  terms  of  £ and  Bx 
(all  of  them  should  be  one  in  an  ideal  situation):  The 
results  are  listed  for  all  three  models  in  column  2 of 
tables  four  to  six. 

Because  of  the  perfect  structure  of  Bj  for  model  2£ 
the  estimates  in  column  2 are  equal  To  the  input  coef- 
ficients. For  method  (1)  and  for  model  2h  the  coef- 
ficients are  reproduced  with  small  distortions  because 
of  the  damped  terms  on  the  diagonal  of  Ay  and  Bj 
respectively.  “ “ 

(2)  The  part  produced  by  the  off-diagonal  terms  of  Ax  and 
Bx  (all  of  them  should  be  zero) : The  results  ar^ given 
in  column  3 of  tables  four  to  six. 

Also  in  column  3 there  are  only  very  small  distortions 
always  less  than  1/1000  of  the  corresponding  coefficient 
for  model  2 1.  The  distrubing  contribution  caused  by  the 
off-diagonal  terms  of  Ax  and  Bx  is  however  significant 
for  method  (1)  and  for* model  Si. 
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(3)  The  part  produced  by  the  matrices  Aa  and  Bg 
which  is  die  alaislng  effect  of  the  coefficienfi  with 
l>  l o on  a,  (all  elements  should  be  zero  in  an  Ideal 
situation).  The  results  are  given  in  column  4 of 
tables  four  to  six. 

The  dimension  of  the  matrices  A,  and  Bs  is  only 
n*(673-k  + l)  instead  of  n*  (® -k  + 1).  'fLe  number 
673  is  thereby  the  full  number  of  potential  coeffi- 
cients up  to  L = 25  and  |m|=  25  excluding  degree 
zero  and  one. 

We  see  in  column  4 that  the  alaising  effect  for  me- 
thod (2)  model  21  is  higher  as  compared  to  method 
(1)  and  also  to  method  (2)  model  2h.  But  since  the 
magnitude  of  the  input  coefficients  decreases  rapidly 
with  increasing  degree,  e.  g.  expressed  by  Kaula's 
rule  of  thumb,  the  disturbing  influence  on  the  esti- 
mated coefficients  is  comparably  small. 


Column  5 of  the  tables  four  to  six  contains  the  sum  of  all  three  parts  which  is  the 
vector  of  the  estimated  coefficients,  a>  and  in  column  6 the  vector,  a^ , of  the 
input  coefficients  is  listed.  The  "best"  estimates  for  the  spherical  harmonic 
coefficients  in  the  sense  that  the  estimated  coefficients  are  closest  to  the  input 
coefficients  is  achieved  by  method  (2)  with  the  low  noise  model  21.  The  reasons 
for  this  fact  shall  be  more  closely  analyzed  in  the  sequel. 

For  convenience  we  denote 
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model  2*« 


0.01  for  method  (3)  modW 


from  off- 
diagonal  of  A, 

from 

5> 

estimate 
sum  2,3.  4 

input 

1 

2 

3 

4 

5 

6 

2 

0 

0.06317 

0 . 166  31 

0.07426 

0.30776 

0.06713 

2 

1 

-0.00041 

-O.OfOll 

-0.06661 

-0.14714 

-0.00042 

2 

1 

0.00180 

-0.07056 

0.03697 

-0.02979 

0.00183 

2 

2 

15.26636 

0.04255 

0.01277 

15.32186 

15.54655 

2 

2 

-0.75019 

o.oeisa 

-0.02123 

-8.68984 

-8.91066 

□ 

5.79696 

0.00009 

•0.01360 

5.78325 

6.11637 

B 

12.70010 

-0.00003 

0.09497 

12.79504 

12.96395 

H 

1.56016 

0.00005 

-0.05432 

1.50  569 

1.59258 

3 

2 

5.60725 

-0.00007 

-0.02939 

5.57779 

5.73012 

3 

2 

-3.09300 

-0.00003 

0.05176 

-3.84126 

-3.97831 

3 

3 

4.44619 

0.00001 

0.00032 

4.44652 

4.57315 

3 

3 

0.79654 

0.00001 

-0.01617 

8.78038 

9.04774 

4 

0 

-1.10642 

-0.00139 

-0.00610 

-1.11391 

-1.59498 

4 

1 

-3.10209 

0.00004 

-0.07775 

-3. 18060 

-3.42616 

4 

1 

-2.73573 

0.00003 

0.04540 

-2.69031 

-3.02075 

4 

2 

2.00410 

0.05235 

0.04531 

2.10184 

2.21423 

4 

2 

3.04163 

-0.02995 

-0.07760 

3.73406 

4.24426 

4 

3 

5.69139 

-0.00004 

0.00032 

5.69168 

6.28301 

4 

3 

-1.13613 

0.00004 

0.04370 

-1.09236 

-1.25423 

4 

4 

-1 .09034 

-0.00004 

-O.OOSOO 

-1.09536 

-1.24628 

4 

4 

1.70339 

0.0 

-0.00690 

1.69448 

1.94700 

Table  4.  Estimated  Geoid  Height  Spherical  Harmonic  Coefficients 
for  Method  (1). 
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from 

diagonal  of  A, 


2 


0.06713 


0.00062 


0.13652 


estimate 
sum  of  2, 3, 6 4 


0.20427 


0.06713 


•0.00042 


0.00046 


-0.02631 


-0.02627 


-0.00042 


0.00183 


0.00048 


0.07752 


0.07982 


0.00183 


15.54655 


-0.00037 


0.05181 


15.59799 


15.54655 


-8.91066 


-0.00069 


0.01823 


-8.89312 


-8.91066 


6.11637 


0.00007 


-0.03422 


6.08221 


6.11637 


12.96395 


-0.00003 


0.14708 


13.11100 


12.96395 


1.59258 


0.00004 


-0.02259 


1.57003 


1.59258 


5.73012 


-3.97831 


4.57315 


-0.00007 


-0.00004 


0.00001 


-0.00068 


0.10914 


0.03952 


5.72937 


-3.86921 


4.61268 


5.73012 


-3.97831 


4.57315 


9.04774 


0.00001 


0.01678 


9.06453 


-1.59498 


-0.00015 


0.06859 


-1.52654 


-1.59498 


-3.42616 


0.00005 


-0.04637 


-3.47249 


-3.42616 


-3.02075 


0.00002 


0.08910 


-2.93163 


-3.02075 


2.21423 


4.24428 


6.2e301 


-0 .00447 


0.00262 


-0 .00005 


0.09618 


-0.05863 


0.03952 


2.30595 


4.18827 


6.32248 


2.21423 


4.24428 


6.28301 


-1.25423 


0.00005 


0.09197 


-1.16222 


-1.25423 


-1.24628 


-0.00005 


0.03340 


-1.21293 


-1.24628 


1.94700 


0.02892 


1.97592 


Table  5.  Estimated  Geoid  Height  Spherical  Harmonic  Coefficients 
for  Method  (2)  Model  Si. 
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from 

diagonal  of  A, 

from  off- 
diagonal  of  A, 

from 

B. 

estimate 
sum  of  2,3,  ft  4 

input 

1 

2 

3 

4 

5 

6 

2 

0 

0.04453 

0.04917 

0.09261 

0.20631 

0.06713 

2 

1 

>0.00040 

-0.14080 

-0.06678 

-0.21599 

-0.00042 

2 

I 

0.00176 

-0.13113 

0.03895 

-0 .09042 

0.00183 

2 

2 

14.79040 

0.05210 

0.01329 

14.86379 

15.54655 

2 

2 

-0.40185 

0.09988 

-0.02234 

-8.40431 

-8.91066 

3 

0 

5.72954 

0.00006 

-0.06883 

5.66070 

6.11637 

3 

1 

12.41922 

-0 .00003 

0.10329 

12.52240 

12.96395 

3 

1 

1.52566 

0.00004 

-0.05925 

1.46645 

1.59258 

3 

2 

5.42104 

-0.00007 

-0.03770 

5.36319 

5.73012 

3 

2 

-3.76372 

-0.00003 

0.06611 

-3.69764 

-3.97831 

3 

3 

4.28032 

0 .00001 

0.00025 

4.28050 

4.57315 

3 

3 

8.46838 

0.00001 

-0.02103 

6.44736 

9.04774 

4 

0 

-1.12950 

-0.00049 

0.02020 

-1.10979 

-1.59498 

4 

1 

-2.92038 

0.00004 

-0.07350 

-2.99364 

-3.42616 

4 

1 

-2.57481 

0.00003 

0.04262 

-2.53196 

-3.02075 

4 

2 

1.83262 

0.06530 

0.04720 

1.94512 

2.21423 

4 

2 

3.51279 

-0.03730 

-0.08108 

3.39433 

4.24428 

4 

3 

5.05440 

-0.00004 

0.00025 

5.05461 

6.28301 

4 

3 

-1.00897 

0.00004 

0.04244 

-0.96649 

-1.25423 

4 

4 

-0.97537 

-0.00004 

-0.00454 

-0.97995 

-1.24628 

4 

4 

1.52378 

0.0 

-0.00804 

1.51574 

1.94700  j 

Table  6.  Estimated  Geoid  Height  Spherical  Harmonic  Coefficients 
for  Method  (2)  Model  2h. 
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where  "A"  indicates  that  the  elements  of  £C  are  usually  smaller  than  these  of 
C,  especially  when  the  maximum  degree  Iq  is  high.  Under  the  usual  conditions 
of  coovergency  we  expand  C^1  = (C  + ACf1  into  a power  series, 

(47)  c;,1  = (C  + AC  I-1  = C'1-  g-xAC  C"1  + g-1  AC  C'1  - ... 

Denoting  the  geoid  height  part  generated  by  the  spherical  harmonic  expansion 
up  to  lo  with  Ui  and  the  residual  geoid  height  with  ua , equation  (31)  becomes 


(48) 


u = R a>  + Sap  = Uj  + up 


Equation  (37)  becomes  with  (47)  and  (48) 

(49)  Sj  - C.u  C Uj  + C.u  C JJji  ■ C,u  C AC  5 1 Uj  - CauC  1 AC  C 1ua  +. . • 
or 

(50)  a>  = C.UC_1  Raj  + 2i«C_1S  Sj,  - C.u  CT1  AC  C_1R  a>  - CauC_1  Ac  C_1S  as  +... 

Multiplying  equation  (50)  on  both  sides  with  • R and  using  equations  (41)  and  (42) 
we  obtain 

(51)  R = R a,  + Sag  - AC  C_1R  a>  - AC  C_1S  a^  + .... 

With  this  series  expansion  it  is  possible  to  discuss  equations  (49)  to  (51) 
term  by  term. 

case  1:  Truncation  of  series  (49)  after  the  first  right  hand 
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term:  Then  the  equation  represents  an  "exact"  collocation 
solution  in  a k -dimensional  space,  i.e.  in  a situation  where 
the  anomalous  potential  of  the  earth  can  be  represented  by 
a series  of  spherical  harmonics  up  to  degree  l0  with  k 
terms.  Equation  (51)  shows  that  truncation  after  the  first 
right  hand  term  results  in  perfect  reproduction  of  a^ . 

case  2:  Truncation  of  series  (49)  to  (51)  after  the  first  two 
terms:  The  first  two  right  hand  terms  of  (50)  are  identical 
to  these  of  equation  (46)  for  method  (2)  if  CBB  = O which 
is  approximately  the  case  for  model  21. 

That  means  that  the  estimation  equations  for  model  21 
with  only  the  first  right  hand  term  reflect  the  situation  des- 
cribed in  case  1,  as  shown  in  table  2 where 


CBU  C"1  R = Bj 


is  a perfect  identity  matrix.  The  second  term  expresses 
the  alaising  effect.  As  is  obvious  from  equations  (50)  and 
(51)  the  magnitude  of  the  elements  of  a^  directly  influences 
the  disturbing  alaising  effect.  Therefore  since  the  magni- 
tude of  the  coefficients  of  quantities  related  to  the  earth's 
anomalous  potential  decreases  in  the  mean  with  increasing 
frequency,  the  alaising  effect  is  comparably  small  as  long 
as  low  harmonic  coefficients  are  estimated. 

case  3:  The  complete  formulas  (49)  to  (51)  represent  the 
collocation  solution  foi  method  (ly.  The  first  and  higher 
order  terms  in  (49)  containing  Uj  and  in  (50)  containing 
aj  cause  the  damping  of  the  diagonal  terms  and  the  appear- 
ance of  off-diagonal  terms  in  matrix  Ax  (table  1).  On  the 
other  hand,  the  first  and  higher  order" terms  containing  Uj, 
and  a e respectively  reduce  significantly  the  alaising  effect 
as  can  be  seen  by  comparing  columns  4 of  tables  4 and  5. 

The  matrix  product  AC  C-1  which  occurs  in  all  first  and 
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higher  order  terms  of  equation  (49)  to  (51)  would  become 
zero  if  instead  of  a finite  point  grid  a global  geoid  height 
function  would  be  given,  for  then  AC  and  C_1  would  be 
orthogonal.  For  this  case  the  solutions  for  method  (1)  and 
method  (2)  would  become  identical  as  already  shown  in 
part  one.  For  all  finite  point  grids  this  orthogonality  rela- 
tionship is  lost  and  causes  the  differences  in  the  solutions 
for  both  methods. 


Using  method  (1)  the  approximation  to  a gravity  function  is  carried  out  in  an 
infinite  dimensional  Hilbert  space  with  kernel  function  C represented  by  an 
infinite  dimensional  Legendre  series.  Therefore  it  would  be  advisable  to 
solve  for  all  coefficients  up  to  infinity.  Only  then  a perfect  reproduction  of 
the  given  geoid  height  from  the  estimated  coefficients  is  ensured. 

Finally,  we  see  that  method  (2)  using  model  2h  yields  almost  the 
same  results  as  method  (1).  Small  differences  may  be  explained  by  the  fact  that 
for  method  (1)  Cuu  is  used  which  contains  all  wavelength  up  to  degree  and  or- 
der infinity,  whereas  for  method  (2)  model  2h  in  equation  (43)  the  second 
right  had  term  i>uu  is  approximated  by 

CO 

AC  = C„  = o*I  with  o3  = 7 s*+1  dt 

Because  of  the  finite  point  grid  the  coefficients  with  l > £0  cannot 
be  resolved.  That  means  that  the  low  wavelength  part  of  the  geoid  heights 
basically  behaves  like  white  noise  or  in  other  words  AC  becomes  approxi- 
mately CBB. 

From  coiumn  5 of  tables  4 to  5 we  see  that  all  three  models  pro- 
duce a reasonable  recovery  of  the  input  geoid  height  spherical  harmonic  coef- 
ficients. Although  method  (2)  with  model  21  produces  the  most  satisfactory 
results  the  underlying  model  can  hardly  be  justified  since  the  low  wavelengths 
components  of  the  given  geoid  heights  are  not  properly  modeled. 
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Part  ID:  Conclusions 

The  two  least  squares  models,  described  by  the  equtions  (1),  (2),  (7), 
and  (8)  can  be  distinguished  by  their  different  variance-covariance  matrices 
Q and  Q*,  compare  equations  (3)  and  (9).  The  differences  in  Q and  Q* 
lead  to  different  formulas  for  the  estimation  of  the  parameters  x,  equations 
(4)  and  (10),  and  for  the  random  signal  s,  equations  (5)  and  (11).  Only  if 
the  relation 

(15)  s'  = Rs 


between  s and  s'  is  valid,  the  estimation  equations  for  x and  s become 
identical.  In  physical  geodesy  where  s'  and  s are  quantities  related  to  the 
earth  anomalous  potential,  equation  (15)  is  only  fulfilled  by  one  of  the  well 
known  operator  equations, such  as  the  Poisson,  Stokes,  or  Vening-Meinesz 
equations.  In  the  frequency  domain,  where  R expresses  a spherical  har- 
monic development  s becomes  the  infinite  dimensional  vector  of  spherical 
harmonic  coefficients  of  this  development.  In  the  parameter  domain,  equa- 
tion (15)  will  be  one  of  the  integral  equations  of  physical  geodesy  and 
we  have  to  assume  s to  be  given  continuously  at  the  entire  sphere,  a. 

For  any  finite  dimension  of  the  random  signal  vector  s we  obtain  dif- 
ferent results  for  x and  s from  Model  One  as  compared  to  Model  Two.  The 
only  approximation  necessary  to  solve  the  finite  problem  is  performed  in  Model 
One  for  the  derivation  of  a continuous  function  interpolated  from  the  observa- 
tions. After  inserting  this  interpolated  function  into  the  operator  equation 
suited  for  the  derivation  of  the  desired  gravity  quantity,  the  operator  equation 
is  solved  without  any  approximation,  whereas  for  Model  Two  the  operator 
equation  itself  is  approximated  by  a finite  matrix  system.  Within  this  respect, 
Model  One  is  superior  to  Model  Two.  The  advantage  of  Model  Two  lies  in  the 
fact,  that  — as  already  mentioned  in  the  introduction  — the  solution  of  a large 
linear  system  with  dimension  equal  to  the  number  of  observations  can  be  re- 
placed by  the  solution  of  a linear  system  with  dimension  equal  to  the  number 
of  unknowns,  since  R appears  explicitly  in  the  solution  equations  (10),  (11) 

) and  allows  the  application  of  certain  matrix  manipulations,  compare  for  example 

(U.  Uotila,  unpublished  paper,  (1973)  or  P.  B.  Liebelt  (1967)  pp. 27). 
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Finally,  any  problem  which  can  be  expressed  by  Model  Two,  i.  e.  for 
which  an  explicit  form  of  the  coefficient  matrix  R can  be  established,  can  also 
be  solved  by  Model  One.  This  is  obvious  since  with  C„  — assumed  to  be  given 
in  both  models — and  R the  covariances  C„»  and  C,V  necessary  in  Model  One 
can  be  derived  from  equations  (13)  and  (14).  The  estimated  parameter  vector 
x and  random  signal  vector  s will  be  identical. 

The  estimated  value  of  a certain  gravity  parameter  depends  not  only 
on  the  special  features  of  the  two  models  but  to  a high  extent  on  the  characteris- 
tics of  the  gravity  field  as  well.  As  an  example  low  degree  and  order  (£  = |m|  = 
4)  spherical  harmonic  coefficients  of  the  geoid  height  have  been  estimated  from 
60  geoid  heights  which  have  been  generated  from  the  GEM-8  set  of  potential 
coefficients  (l  - |m[  =25)  and  were  given  in  a global  regular  grid  over  the 
sphere.  For  this  problem  the  estimation  equations  become  almost  identical 
for  both  methods  with  the  only  difference  that  for  Method  One  the  covariances 
between  the  observations  contain  all  wavelengths  inherent  also  in  the  spectrum 
of  the  observations  whereas  these  of  Method  Two  contain  only  the  wavelength 
range  of  the  coefficients  to  be  estimated.  Because  of  this  a considerable  a- 
laising  effect  would  be  expected  for  Method  Two.  But  the  results  show  that 
because  of  the  rapid  decrease  of  the  degree  variances  of  geoid  heights  with  in- 
creasing degree  the  disturbing  effect  of  the  shorter  wavelengths  (4  -ft*  25)  is 
almost  negligible.  The  best  recovery  of  the  harmonic  coefficients  was  achieved 
using  Model  Two  and  assuming  that  the  given  geoid  heights  contain  no  wavelengths 
with  l> 4.  Method  Two  modeling  the  wavelengths  between  l = 5 and  l = 25  as 
white  noise  yielded  about  the  same  results  as  Method  One. 
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